home *** CD-ROM | disk | FTP | other *** search
/ Computer Select (Limited Edition) / Computer Select.iso / dobbs / v17n06 / contour.asc < prev    next >
Encoding:
Text File  |  1992-05-18  |  18.8 KB  |  694 lines
  1. _CONTOURING DATA FIELDS_
  2. by Bruce (Bear) Giles
  3.  
  4. [LISTING ONE]
  5.  
  6. #include    <stdio.h>
  7. #include    <math.h>
  8. #include    <malloc.h>
  9. #include    <values.h>
  10.  
  11. #if defined (NEVER)
  12. #include    <ieeefp.h>
  13. #else
  14. #define     NaN         0xFFFFFFFF
  15. #define     isnanf(x)   ((x) == NaN)
  16. #endif
  17.  
  18. typedef unsigned    short   ushort;
  19. typedef unsigned    char    uchar;
  20.  
  21. #define DEFAULT_LEVELS  16
  22.  
  23. /* Mnemonics for contour line bearings  */
  24. #define EAST        0
  25. #define NORTH       1
  26. #define WEST        2
  27. #define SOUTH       3
  28.  
  29. /* Mnemonics for relative data point positions  */
  30. #define SAME        0
  31. #define NEXT        1
  32. #define OPPOSITE    2
  33. #define ADJACENT    3
  34.  
  35. /* Bit-mapped information in 'map' field.  */
  36. #define EW_MAP      0x01
  37. #define NS_MAP      0x02
  38.  
  39. typedef struct
  40.     {
  41.     float   x, y;
  42.     }   LIST;
  43. typedef struct
  44.     {
  45.     short   dim_x;          /*  dimensions of grid array... */
  46.     short   dim_y;
  47.     float   max_value;      /*  statistics on the data...   */
  48.     float   min_value;
  49.     float   mean;
  50.     float   std;
  51.     short   contour_mode;   /*  control variable            */
  52.     float   first_level;    /*  first (and subsequent)      */
  53.     float   step;           /*    contour level             */
  54.     char    format[20];     /*  format of contour labels    */
  55.     float   *data;          /*  pointer to grid data        */
  56.     char    *map;           /*  pointer to "in-use" map     */
  57.     LIST    *list;          /*  used by 'Polyline()'        */
  58.     ushort  count;
  59.     }   GRID;
  60. typedef struct
  61.     {
  62.     short   x;
  63.     short   y;
  64.     uchar   bearing;
  65.     }   POINT;
  66. #define MXY_to_L(g,x,y)     ((ushort) (y) * (g)->dim_x + (ushort) (x) + 1)
  67. #define XY_to_L(g,x,y)      ((ushort) (y) * (g)->dim_x + (ushort) (x))
  68.  
  69. extern  void    Text();
  70. extern  void    Polyline();
  71.  
  72. /* Contour generation. */
  73. void    Contour();
  74.  
  75. int     scaleData();
  76. void    startLine();
  77.  
  78. void    startInterior();
  79. void    drawLine();
  80.  
  81. void    markInUse();
  82. uchar   faceInUse();
  83. float   getDataPoint();
  84.  
  85. void    initPoint();
  86. void    lastPoint();
  87. uchar   savePoint();
  88.  
  89. /* inc     > 0 increment to use between contour levels.
  90.  *         < 0 number of contour levels to generate [abs(inc)].
  91.  *         = 0 generate default number of contour levels.
  92.  */
  93. void    Contour (data, dim_x, dim_y, inc)
  94. float   *data;
  95. int     dim_x;
  96. int     dim_y;
  97. double  inc;
  98.     {
  99.     GRID    grid;
  100.  
  101.     grid.data   = data;
  102.     grid.dim_x  = dim_x;
  103.     grid.dim_y  = dim_y;
  104.  
  105.     /* Allocate buffers used to contain contour information */
  106.     if ((grid.map = malloc ((dim_x + 1) * dim_y)) == NULL)
  107.         {
  108.         fprintf (stderr, "Contour(): unable to allocate buffer! (%d bytes)\n",
  109.             (dim_x + 1) * dim_y * sizeof (LIST));
  110.  
  111.         free ((char *) grid.map);
  112.         return;
  113.         }
  114.     grid.list = (LIST *) malloc (2 * dim_x * dim_y * sizeof (LIST));
  115.     if (grid.list == (LIST  *) NULL)
  116.         {
  117.         fprintf (stderr, "Contour(): unable to allocate buffer! (%d bytes)\n",
  118.             2 * dim_x * dim_y * sizeof (LIST));
  119.  
  120.         free ((char *) grid.map);
  121.         return;
  122.         }
  123.     /* Generate contours, if not a uniform field. */
  124.     if (scaleData (&grid, inc))
  125.  
  126.         startLine (&grid);
  127.     /* Release data structures. */
  128.     free ((char *) grid.map);
  129.     free ((char *) grid.list);
  130.     }
  131.  
  132. /* scaleData--Determine necessary statistics for contouring data set: global 
  133.  * maximum & minimum, etc. Then initialize items used by rest of module. */
  134. int     scaleData (grid, inc)
  135. GRID    *grid;
  136. double  inc;
  137.     {
  138.     ushort  i;
  139.     float   step, level;
  140.     float   sum, sum2, count;
  141.     float   p, *u, *v, r;
  142.     char    *s;
  143.     short   n1, n2;
  144.     int     first, n;
  145.     long    x;
  146.  
  147.     sum = sum2 = count = 0.0;
  148.  
  149.     first = 1;
  150.     s = grid->map;
  151.     u = grid->data;
  152.     v = u + grid->dim_x * grid->dim_y;
  153.     for (i = 0; i < grid->dim_x * grid->dim_y; i++, u++, v++, s++)
  154.         {
  155.         r = *u;
  156.         sum   += r;
  157.         sum2  += r * r;
  158.         count += 1.0;
  159.  
  160.         if (first)
  161.             {
  162.             grid->max_value = grid->min_value = r;
  163.             first = 0;
  164.             }
  165.         else if (grid->max_value < r)
  166.             grid->max_value = r;
  167.         else if (grid->min_value > r)
  168.             grid->min_value = r;
  169.         }
  170.     grid->mean = sum / count;
  171.     if (grid->min_value == grid->max_value)
  172.         return 0;
  173.     grid->std  = sqrt ((sum2 - sum * sum /count) / (count - 1.0));
  174.     if (inc > 0.0)
  175.         {
  176.         /* Use specified increment */
  177.         step = inc;
  178.         n = (int) (grid->max_value - grid->min_value) / step + 1;
  179.  
  180.         while (n > 40)
  181.             {
  182.             step *= 2.0;
  183.             n = (int) (grid->max_value - grid->min_value) / step + 1;
  184.             }
  185.         }
  186.     else 
  187.         {
  188.         /* Choose [specified|reasonable] number of levels and normalize 
  189.          * increment to a reasonable value. */
  190.         n = (inc == 0.0) ? DEFAULT_LEVELS : (short) fabs (inc);
  191.  
  192.         step = 4.0 * grid->std / (float) n;
  193.         p    = pow (10.0, floor (log10 ((double) step)));
  194.         step = p * floor ((step + p / 2.0) / p);
  195.         }
  196.     n1 = (int) floor (log10 (fabs (grid->max_value)));
  197.     n2 = -((int) floor (log10 (step)));
  198.  
  199.     if (n2 > 0)
  200.         sprintf (grid->format, "%%%d.%df", n1 + n2 + 2, n2);
  201.     else 
  202.         sprintf (grid->format, "%%%d.0f", n1 + 1);
  203.     if (grid->max_value * grid->min_value < 0.0)
  204.         level = step * floor (grid->mean / step);           /*  odd */
  205.     else
  206.         level = step * floor (grid->min_value / step);
  207.     level -= step * floor ((float) (n - 1)/ 2);
  208.  
  209.     /* Back up to include add'l levels, if necessary */
  210.     while (level - step > grid->min_value)
  211.         level -= step;
  212.  
  213.     grid->first_level = level;
  214.     grid->step = step;
  215.     return  1;
  216.     }
  217. /* startLine -- Locate first point of contour lines by checking edges of 
  218.    gridded data set, then interior points, for each contour level. */
  219. static
  220. void    startLine (grid)
  221. GRID    *grid;
  222.     {
  223.     ushort  idx, i, edge;
  224.     double  level;
  225.     for (idx = 0, level = grid->first_level; level < grid->max_value; 
  226.         level += grid->step, idx++)
  227.         {
  228.  
  229.         /* Clear flags */
  230.         grid->contour_mode = (level >= grid->mean);
  231.         memset (grid->map, 0, grid->dim_x * grid->dim_y);
  232.  
  233.         /* Check edges */
  234.         for (edge = 0; edge < 4; edge++)
  235.             startEdge (grid, level, edge);
  236.         /* Check interior points */
  237.         startInterior (grid, level);
  238.         }
  239.     }
  240. /* startEdge -- For a specified contour level and edge of gridded data set,
  241.  *  check for (properly directed) contour line. */
  242. static
  243. void    startEdge (grid, level, bearing)
  244. GRID    *grid;
  245. float   level;
  246.  
  247. uchar   bearing;
  248.     {
  249.     POINT   point1, point2;
  250.     float   last, next;
  251.     short   i, ds;
  252.  
  253.     switch (point1.bearing = bearing)
  254.         {
  255.         case EAST:
  256.             point1.x = 0;
  257.             point1.y = 0;
  258.             ds = 1;
  259.             break;
  260.         case NORTH:
  261.             point1.x = 0;
  262.             point1.y = grid->dim_y - 2;
  263.             ds = 1;
  264.             break;
  265.         case WEST:
  266.             point1.x = grid->dim_x - 2;
  267.             point1.y = grid->dim_y - 2;
  268.             ds = -1;
  269.             break;
  270.         case SOUTH:
  271.             point1.x = grid->dim_x - 2;
  272.             point1.y = 0;
  273.             ds = -1;
  274.             break;
  275.         }
  276.     switch (point1.bearing)
  277.         {
  278.             /* Find first point with valid data. */
  279.         case EAST:
  280.         case WEST:
  281.             next = getDataPoint (grid, &point1, SAME);
  282.             memcpy ((char *) &point2, (char *) &point1, sizeof (POINT));
  283.             point2.x -= ds;
  284.  
  285.             for (i = 1; i < grid->dim_y; i++, point1.y = point2.y += ds)
  286.                 {
  287.                 last = next;
  288.                 next = getDataPoint (grid, &point1, NEXT);
  289.                 if (last >= level && level > next)
  290.                    {
  291.                    drawLine (grid, &point1, level);
  292.                    memcpy ((char *) &point1, (char *) &point2, sizeof (POINT));
  293.                    point1.x = point2.x + ds;
  294.  
  295.                    }
  296.                 }
  297.             break;
  298.             /* Find first point with valid data. */
  299.         case SOUTH:
  300.         case NORTH:
  301.             next = getDataPoint (grid, &point1, SAME);
  302.             memcpy ((char *) &point2, (char *) &point1, sizeof (POINT));
  303.             point2.y += ds;
  304.  
  305.             for (i = 1; i < grid->dim_x; i++, point1.x = point2.x += ds)
  306.                 {
  307.                 last = next;
  308.                 next = getDataPoint (grid, &point1, NEXT);
  309.                 if (last >= level && level > next)
  310.                     {
  311.                    drawLine (grid, &point1, level);
  312.  
  313.                    memcpy ((char *) &point1, (char *) &point2, sizeof (POINT));
  314.                    point1.y = point2.y - ds;
  315.                    }
  316.                 }
  317.             break;
  318.         }
  319.     }
  320. /* startInterior -- For a specified contour level, check for (properly 
  321. directed) contour line for all interior data points. Do _not_ follow contour 
  322. lines detected by the 'startEdge' routine. */
  323. static
  324. void    startInterior (grid, level)
  325. GRID    *grid;
  326. float   level;
  327.     {
  328.     POINT   point;
  329.     ushort  x, y;
  330.     float   next, last;
  331.     for (x = 1; x < grid->dim_x - 1; x++)
  332.         {
  333.         point.x = x;
  334.         point.y = 0;
  335.         point.bearing = EAST;
  336.         next = getDataPoint (grid, &point, SAME);
  337.         for (y = point.y; y < grid->dim_y; y++, point.y++)
  338.  
  339.             {
  340.             last = next;
  341.             next = getDataPoint (grid, &point, NEXT);
  342.             if (last >= level && level > next)
  343.                 {
  344.                 if (!faceInUse (grid, &point, WEST))
  345.                     {
  346.                     drawLine (grid, &point, level);
  347.                     point.x = x;
  348.                     point.y = y;
  349.                     point.bearing = EAST;
  350.                     }
  351.                 }
  352.             }
  353.         }
  354.     }
  355. /* drawLine -- Given in initial contour point by either 'startEdge' or 
  356. 'startInterior', follow contour line until it encounters either an edge 
  357. or previously contoured cell. */
  358. static
  359. void    drawLine (grid, point, level)
  360. GRID    *grid;
  361. POINT   *point;
  362. float   level;
  363.     {
  364.     uchar   exit_bearing;
  365.     uchar   adj, opp;
  366.     float   fadj, fopp;
  367.  
  368.     initPoint (grid);
  369.  
  370.     for ( ;; )
  371.         {
  372.         /* Add current point to vector list. If either of the points is 
  373.          * missing, return immediately (open contour). */
  374.         if (!savePoint (grid, point, level))
  375.             {
  376.             lastPoint (grid);
  377.             return;
  378.             }
  379.         /* Has this face of this cell been marked in use? If so, then this is 
  380.          * a closed contour. */
  381.         if (faceInUse (grid, point, WEST))
  382.             {
  383.             lastPoint (grid);
  384.             return;
  385.             }
  386.         /* Examine adjacent and opposite corners of cell; determine 
  387.          * appropriate action. */
  388.         markInUse (grid, point, WEST);
  389.  
  390.         fadj = getDataPoint (grid, point, ADJACENT);
  391.         fopp = getDataPoint (grid, point, OPPOSITE);
  392.  
  393.         /* If either point is missing, return immediately (open contour). */
  394.         if (isnanf (fadj) || isnanf (fopp))
  395.             {
  396.             lastPoint (grid);
  397.             return;
  398.             }
  399.         adj = (fadj <= level) ? 2 : 0;
  400.         opp = (fopp >= level) ? 1 : 0;
  401.         switch (adj + opp)
  402.             {
  403.             /* Exit EAST face. */
  404.             case 0: 
  405.                 markInUse (grid, point, NORTH);
  406.                 markInUse (grid, point, SOUTH);
  407.                 exit_bearing = EAST;
  408.                 break;
  409.             /* Exit SOUTH face. */
  410.             case 1:
  411.                 markInUse (grid, point, NORTH);
  412.                 markInUse (grid, point, EAST);
  413.                 exit_bearing = SOUTH;
  414.                 break;
  415.             /* Exit NORTH face. */
  416.             case 2:
  417.                 markInUse (grid, point, EAST);
  418.  
  419.                 markInUse (grid, point, SOUTH);
  420.                 exit_bearing = NORTH;
  421.                 break;
  422.             /* Exit NORTH or SOUTH face, depending upon contour level. */
  423.             case 3:
  424.                 exit_bearing = (grid->contour_mode) ? NORTH : SOUTH;
  425.                 break;
  426.             }
  427.         /* Update face number, coordinate of defining corner. */
  428.         point->bearing = (point->bearing + exit_bearing) % 4;
  429.         switch (point->bearing)
  430.             {
  431.             case EAST : point->x++; break;
  432.             case NORTH: point->y--; break;
  433.             case WEST : point->x--; break;
  434.             case SOUTH: point->y++; break;
  435.             }
  436.         }
  437.     }
  438. /* markInUse -- Mark the specified cell face as contoured. This is necessary to
  439.  *  prevent infinite processing of closed contours. see also:  faceInUse */
  440. static
  441. void    markInUse (grid, point, face)
  442. GRID    *grid;
  443. POINT   *point;
  444. uchar   face;
  445.     {
  446.     face = (point->bearing + face) % 4;
  447.     switch (face)
  448.         {
  449.         case NORTH:
  450.         case SOUTH:
  451.             grid->map[MXY_to_L (grid,
  452.                 point->x, point->y + (face == SOUTH ? 1 : 0))] |= NS_MAP;
  453.             break;
  454.         case EAST:
  455.         case WEST:
  456.             grid->map[MXY_to_L (grid,
  457.                 point->x + (face == EAST ? 1 : 0), point->y)] |= EW_MAP;
  458.  
  459.             break;
  460.         }
  461.     }
  462. /* faceInUse -- Determine if the specified cell face has been marked as 
  463.  * contoured. This is necessary to prevent infinite processing of closed 
  464.  * contours. see also: markInUse */
  465. static
  466. uchar   faceInUse (grid, point, face)
  467. GRID    *grid;
  468. POINT   *point;
  469. uchar   face;
  470.     {
  471.     uchar   r;
  472.     face = (point->bearing + face) % 4;
  473.     switch (face)
  474.         {
  475.         case NORTH:
  476.         case SOUTH:
  477.             r = grid->map[MXY_to_L (grid,
  478.                     point->x, point->y + (face == SOUTH ? 1 : 0))] & NS_MAP;
  479.             break;
  480.         case EAST:
  481.         case WEST:
  482.             r = grid->map[MXY_to_L (grid,
  483.                     point->x + (face == EAST ? 1 : 0), point->y)] & EW_MAP;
  484.             break;
  485.         }
  486.     return r;
  487.     }
  488.  
  489.  
  490. /* initPoint -- Initialize the contour point list. 
  491.  * see also: savePoint, lastPoint  */
  492. static
  493. void    initPoint (grid)
  494. GRID    *grid;
  495.     {
  496.     grid->count = 0;
  497.     }
  498.  
  499. /* lastPoint -- Generate the actual contour line from the contour point list.
  500.  *  see also:  savePoint, initPoint  */
  501. static
  502. void    lastPoint (grid)
  503. GRID    *grid;
  504.     {
  505.     if (grid->count)
  506.         Polyline (grid->count, grid->list);
  507.     }
  508.  
  509. /* savePoints -- Add specified point to the contour point list.
  510.  *  see also:  initPoint, lastPoint */
  511. static  
  512. uchar   savePoint (grid, point, level)
  513. GRID    *grid;
  514. POINT   *point;
  515. float   level;
  516.     {
  517.     float   last, next;
  518.     float   x, y, ds;
  519.     char    s[80];
  520.  
  521.     static  int cnt = 0;
  522.  
  523.     last = getDataPoint (grid, point, SAME);
  524.     next = getDataPoint (grid, point, NEXT);
  525.  
  526.     /* Are the points the same value? */
  527.     if (last == next)
  528.         {
  529.         fprintf (stderr, "(%2d, %2d, %d)  ", point->x, point->y, 
  530.             point->bearing);
  531.         fprintf (stderr, "%8g  %8g  ", last, next);
  532.         fprintf (stderr, "potential divide-by-zero!\n");
  533.         return 0;
  534.         }
  535.     x = (float) point->x;
  536.     y = (float) point->y;
  537.     
  538.     ds = (float) ((last - level) / (last - next));
  539.  
  540.  
  541.     switch (point->bearing)
  542.         {
  543.         case EAST :                 y += ds;        break;
  544.         case NORTH: x += ds;        y += 1.0;       break;
  545.         case WEST : x += 1.0;       y += 1.0 - ds;  break;
  546.         case SOUTH: x += 1.0 - ds;                  break;
  547.         }
  548.  
  549.     /* Update to contour point list */
  550.     grid->list[grid->count].x = x / (float) (grid->dim_x - 1);
  551.     grid->list[grid->count].y = y / (float) (grid->dim_y - 1);
  552.  
  553.     /* Add text label to contour line. */
  554.     if (!(cnt++ % 11))
  555.         {
  556.         sprintf (s, grid->format, level);
  557.         Text (s, grid->list[grid->count].x, grid->list[grid->count].y);
  558.         }
  559.     /* Update counter */
  560.     grid->count++;
  561.  
  562.     return 1;
  563.     }
  564. /* getDataPoint -- Return the value of the data point in specified corner of 
  565.  * specified cell (the 'point' parameter contains the address of the
  566.  *  top-left corner of this cell). */
  567. static  
  568. float   getDataPoint (grid, point, corner)
  569. GRID    *grid;
  570. POINT   *point;
  571. uchar   corner;
  572.     {
  573.     ushort  dx, dy;
  574.     ushort  offset;
  575.  
  576.     switch ((point->bearing + corner) % 4)
  577.         {
  578.         case SAME    :  dx = 0; dy = 0; break;
  579.         case NEXT    :  dx = 0; dy = 1; break;
  580.         case OPPOSITE:  dx = 1; dy = 1; break;
  581.         case ADJACENT:  dx = 1; dy = 0; break;
  582.  
  583.         }
  584.     offset = XY_to_L (grid, point->x + dx, point->y + dy);
  585.     if ((short) (point->x + dx) >= grid->dim_x || 
  586.         (short) (point->y + dy) >= grid->dim_y || 
  587.         (short) (point->x + dx) < 0 ||
  588.         (short) (point->y + dy) < 0)
  589.         {
  590.         return  NaN;
  591.         }
  592.     else
  593.         return grid->data[offset];
  594.     }
  595.  
  596.  
  597.  
  598. [LISTING TWO]
  599.  
  600. #include    <stdio.h>
  601.  
  602. typedef struct
  603.     {
  604.     float   x, y;
  605.     }   LIST;
  606. void    Polyline (n, list)
  607. int     n;
  608. LIST    *list;
  609.     {
  610.     short   x0, x1, y0, y1; 
  611.     if (n < 2)
  612.         return;
  613.     printf ("newpath\n");
  614.     printf ("%.6f  %.6f moveto\n", list->x, 1.0 - list->y);
  615.     list++;
  616.  
  617.     while (--n)
  618.         {
  619.         printf ("%.6f  %.6f lineto\n", list->x, 1.0 - list->y);
  620.         list++;
  621.         }
  622.     printf ("stroke\n\n");
  623.  
  624.     }
  625. void    Text (s, x, y)
  626. char    *s;
  627. float   x, y;
  628.     {
  629.     printf ("%.6f  %.6f moveto (%s) show\n", x, 1.0 - y, s);
  630.     }
  631. void    psOpen ()
  632.     {
  633.     printf ("%%!\n");
  634.     printf ("save\n");
  635.     printf ("\n");
  636.     printf ("/Helvetica findfont 0.015 scalefont setfont\n");
  637.     printf ("\n");
  638.  
  639.     printf ("72 252 translate\n");
  640.     printf ("468 468 scale\n");
  641.     printf ("0.001 setlinewidth\n");
  642.     printf ("\n");
  643.     printf ("newpath\n");
  644.     printf (" 0 0 moveto\n");
  645.     printf (" 0 1 lineto\n");
  646.     printf (" 1 1 lineto\n");
  647.     printf (" 1 0 lineto\n");
  648.     printf (" closepath\n");
  649.     printf ("stroke\n");
  650.     printf ("clippath\n");
  651.     printf ("\n");
  652.     printf ("0.00001 setlinewidth\n");
  653.     printf ("\n");
  654.     }
  655. void    psClose ()
  656.     {
  657.     printf ("restore\n");
  658.     printf ("showpage\n");
  659.     }
  660.  
  661.  
  662. [LISTING THREE]
  663.  
  664. #include    <stdio.h>
  665. #include    <math.h>
  666.  
  667. float   data[400];
  668.  
  669. void    main (argc, argv)
  670. int     argc;
  671. char    **argv;
  672.     {
  673.     float   *s;
  674.     int     i, j;
  675.     double  x, y, r1, r2;
  676.  
  677.     s = data;
  678.     for (i = 0, s = data; i < 20; i++)
  679.         for (j = 0; j < 20; j++)
  680.             {
  681.             x = ((double) i - 9.5) / 4.0;
  682.             y = ((double) j - 9.5) / 4.0;
  683.  
  684.             *s++ = (float) (10.0 * cos(x) * cos(y));
  685.             }
  686.     psOpen ();
  687.     Contour (data, 20, 20, 2.0);
  688.  
  689.     psClose ();
  690.     }
  691.  
  692.  
  693.  
  694.